Motion of a random walker in a quenched power law correlated velocity field 
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We study the motion of a random walker in one longitudinal and d transverse dimensions 
with a quenched power law correlated velocity field in the longitudinal x-direction. The model 
is a modification of the Matheron-de Marsily (MdM) model, with long-range velocity correlation. 
For a velocity correlation function, dependent on transverse co-ordinates y as l/(a -|- \yi — y2\)°', 
we analytically calculate the two-time correlation function of the a;-coordinate. We find that 
the motion of the a;-coordinate is a fractional Brownian motion (fBm), with a Hurst exponent 
H — max[l/2, (1 — a/4), (1 — d/A)]. From this and known properties of fBM, we calculate the dis- 
order averaged persistence probability of x{t) up to time t. We also find the lines in the parameter 
space of d and a along which there is marginal behaviour. We present results of simulations which 
support our analytical calculation. 

PACS numbers: 02.50.Ey, 05.40.-a, 05.60.-k 



I. INTRODUCTION 



It is well known that a random walker making short-range hops, namely a Brownian particle in a thermal envi- 
ronment, has an average incremental displacement which goes as square root of the time difference. Although this 
property is rather robust, deviation from this diffusive behaviour has been observed in a wide variety of physical 
situations T*] . Super-diffusion is well known in the context of turbulence and Levy walks |2| ■ Yet another way to see 
anomalous diffusion is to subject the random walker to a quenched velocity field. A concrete example of this is the 
Matheron-de Marsily (MdM) model, which was introduced to study the hydrodynamic dispersion of solute particles 
in sedimentary layered porous rock formation In the MdM model a single particle diffuses in a (1 -f l)-dimensional 
layered medium, such that the disorder is anisotropic. The motion of the particle is purely Brownian along the 
transverse y direction, while along the longitudinal x direction, added to the thermal noise, the particle is driven by 
a quenched drift velocity v{y) that is a random function of only the transverse coordinate y. The coupling between 
diffusion and convection due to the spatially random but temporally static velocity field, generates a typical bias in 
the X direction giving rise to a super-diffusive longitudinal transport for large t P, 0, 0, 01 • 

The temporal behaviour is quantitatively described by the standard two-time dependent correlation function, and 
the persistence probability. Intuitively, persistence of a system concerns the property of the system to stay in some 
given state. The persistence probability P(t) which is related to the first-passage probability F(t), as F{t) = —dP{t)/dt 
is simply the probability that the particle does not cross a given point up to time t. In systems with super-diffusion, 
P{t) is expected to get enhanced. Persistence properties in disordered systems are harder to find, and thus results 
known for such models like the Sinai model 0,11, 1^ and the MdM model P, IToj are of great importance. Recently 
for a {d+1) dimensional version of MdM model the persistence probability was analytically derived in pd| . 

In earlier studies of the MdM model [111 in 113, only short range correlated v{y) had been considered. However, 
since originally the model was motivated [3 by general flow in certain types of fractured rocks with parallel fractures 
which allow the propagation of dissolved species, it is possible that the flows parallel to the bedding might have 
velocity correlations which are spread over long distances. It may be mentioned that different aspects of dynamics in 
MdM fields with algebraic correlations, have been studied Generally long range correlated disordered problems 
are hard to solve. For isotropic disorder with power law correlations, various results are known [ll [T3llT3 |. The MdM 
model being anisotropic allows for exact anatytic treatment as compared to field theoretic treatment done in other 
cases. In this paper we have been able to analytically derive the incremental correlation function and persistence 
probability, for the MdM model with long range power law correlated v[y). We find that long range correlation 
incorporates new temporal behaviour, as compared to the short range MdM model. The relevant exponent of P{t) 
becomes dependent on the power law exponent a and the transverse dimension d. 
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The method that we follow is similar to 0, and was used even earlier in ^E|- We first show that x{t) is 
a fractional Brownian motion (fBm), and then use its known first passage property to solve for the persistence 
probability. A stochastic process x(t) (with zero mean {x{t)) = 0) is called a fBm if its incremental two-time 
correlation function C{ti,t2) = {{x{ti) — x{t2))'^) is stationary, i.e., depends only on the difference \ti — and 
moreover, grows asymptotically as a power law, 

Cih,t2)r^\ti-t2\''", |tl-i2|»l. (1) 

The parameter H is called the Hurst exponent that characterizes the fBm'lT'l and (. . .) denotes the expectation value 
over all realizations of the process. For example, the ordinary Brownian motion is a fBm with H = 1/2. The zero 
crossing properties of fBm has been studied before[la,llMliEl20| and it is known, both analytically and numerically, 
that the probability that a fBm does not cross zero up to time t has a power law decay, P{t) ^ for large t with 
9=1 — H^l&^,2Ql■ We use this property to derive 9 from H. 



II. MODEL AND MAIN RESULTS 



We consider a generalized version of the MdM model in d{> 0) transverse plus 1 longitudinal spatial dimensions. 
The transverse coordinates of the particle, denoted by j/i(t)'s with i = 1, 2, . . . , d, undergo ordinary Brownian motion, 

y^^V^it), (2) 

where i]i's are zero mean Gaussian white noises with correlators, {i]i{'t)Vj{t')) — — t'). Equivalently we can use 

a single equation involving d-dimensional vectors y = {yi{t)} and ff = {r]i{t)}: 

= (3) 

The longitudinal coordinate x{t), in contrast, is driven by a random drift v[y{t)] that depends only on the transverse 
coordinates 

x = v[yit)]+vAt), (4) 

where 77a; is again a zero mean Gaussian white noise with {rjx{ti)7]x{t2)) ~ — ^2)- The velocity field v[y\ is quenched, 
i.e., constant in time. Operationally this means that for a given realization of the function v, one has to evolve x(t) 
over 'thermal histories' in Eq. and then one needs to 'disorder average' (denoted by t^) over different realizations 
of the random function v. In the original MdM model and in the latter works on slightly modified versions of it 
Qi IE ^ IE IE ^1 the random function v{y) has been considered to be a Gaussian with zero mean {v{y) ^ 0), 

with a correlator v(jji)v{y2), which is either S{yi — j/2) or is short ranged namely l/[27rda^]'^/^ exp — (yi — 2/2)^/2^0 

(where a is a small cutoff). 

In contrast to the above, here we consider the random function v to have long-range power-law correlation: 



viy) - 
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v{yi)v{y2) = 7 — -p ^T— . (5) 

[a + m - y2\r 

The exponent a > denotes the strength of the correlation function, and the short distance cutoff scale 'a' ensures 
that there is no divergence as \yi — y2\ ^ 0. This model is analytically treated in Section IlIII and numerically in 
Section Hvl We will see that the role played by the scale a is very interesting. For a < d, it is irrelevant and can be 
set to zero, but for a > d, it has to be finite to ensure a well-behaved asymptotic behaviour. 

Our main result is that the thermal and disorder averaged two-time correlation function of the a:-coordinate, 

{{x{ti) - x{t2)f) ~ \t2 - iiP^. Thus the process x{t) is an fBM. The Hurst exponent H = max[l/2, (1 - a/4), (1 - 
(i/4)]. This means that for d < 2, the exponent will change from iJ = (I — a/4) for a < d, to = (1 — d/4) for 
a > d. On the other hand, for d > 2, iJ = (f — a/4) for a < 2, and H = 1/2 for a > 2. It follows from this that the 
persistence probability of P{t) ~ t^^ , with 9 = 1 — H = min[f /2, a/4, d/4]. We have also found the lines in d and a 
plane along which there is marginal behaviour. 
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III. TWO-TIME CORRELATION FUNCTION AND PERSISTENCE PROBABILITY 

The analytical treatment that we present here is similar to that of |llL Il5j | up to the point that the explicit power 
law form of velocity correlation is substituted. Nevertheless, for completeness, we give the steps briefly. In order to 
solve for correlation functions of x{t) from Eq. ^ we need knowledge of correlation functions of y{t). By integrating 
Eq. 13 we get 



f]{T)dT (6) 

From the above linear relation we conclude that since ^ is a Gaussian random variable, so is y. Moreover another 
vector Y — ?/(ti) — y(T2) which will be useful in the following calculations is also a Gaussian random variable by virtue 
of the linear relation. Using Eq. IHl 



Y = / riiT)dT - / rf{T)dT, (7) 
Jq Jo 

and it follows that (Y) — and variance (y^) = din — T2\, since {'r)i{t)r]j{t')) — Sij6{t — t'). Thus explicitly the 
distribution 

P{Y) = [27rd|ri - T2\]-^/^ exp [-fV2d|ri - T2I] . (8) 

With the above solution of the transverse coordinates y{t), we now proceed to treat the longitudinal coordinate 
x{t). We note that the averages so far over the transverse coordinates, namely (. . .) were thermal averages. The 
coordinate x{t) depends on both a thermal noise rj^ and a quenched random noise v. So we need to do a thermal as 
well as quenched disorder average, which we will denote by (...), where in particular the denotes the quenched 
disorder average. Integrating Eq. we get, 

xit)= I Vx(r)dT+ f v[y{T)]dT. (9) 



From Eq. © we get {x{t)) = 0, while using {•qx{ti)^x{t2)) = 5{ti — t2), we get 



{xiti)xit2)) = min(ti, ia) + /(^i, ^2), (10) 

where the integral I{ti, 12) is given by 



tl /-t2 



I{tl,t2) - / / {v[y{T2)]v[y{T2)])dTldT2 

Jo Jo 
rti rt2 / I \ 

^ )dTldT2. (11) 

In the above, we have used the correlator of v from Eq. Eland Y is the same as in Eq. [3 Since we know the probability 
distribution P{Y) explicitly from Eq. |S| we can proceed from Eq. ^]to evaluate 

tl — 



I(ti,t2) = dTidT2— ■ 

^ ' '0 Jo J-oo {\Y\ + ar{27rd\T, ~ T2\r/^ 

2 rti rt2 rco dTidT2e-^^' Z'^-^dZ 



(Zv/2d|ri -T2I +a) 
e-^'z'^-^ J{Z,ti,t2) dZ. (12) 



r(d/2) Jo 



In Eq. El we have substituted \Y\ = Zy/2d\Ti — T2I in going from the first line to the second, and then performed 
the two time integrals exactly to get explicitly the function 

Jo Jo {Zy/2d\Ti-T2\+a)'- 
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1 2a2-°min(ti,i2) (a + ZV^M)^-" + (a + ZV2M)'^~" - ~ (a + ^^2^^1*2 - ^i| 



(2-a)(l-a) Z4d2(2-a)(4-a) 



+ 



(3 - 2a) -a(a + ZySM)^-" - a{a + Z^/2t^f-" + a^-" + a{a + Z ^2d\t2 - h])^-'^ 



Although Eq. E|has many terms, the two time-correlation fmiction C(fi,<2) = (^^^(ii)) + {x'^it^)) ~ '^{x{ti)x{t2)) 
has a dependence only on t = |ti — t2|j a-nd has the following comparatively simpler expression: 

C(i)=i+jT^^y dZ e-^'z'^-^{^ + ^ [ 2B{a + zV2dt)^-'' ~2Ca{a + zV2dif'" 

+ 2Da^{a + zVmf-" - {2B -2C + 2D)a^-"]} (15) 

In Eq. CHI we have used the symbols A = 2a^-°'/[d(2 - a){l - a)], B ^ 1/M^(2 - a) (4 - a)], C =_(3 - 2a)/[d^{l - 
a){2 - a)(3 - a)] and D = - a)(2 - a)]. Further, Eq. ^can be written as C{t) ^ t + C, where C is the 

integral which has to be analyzed to get the leading temporal behaviour. 

Defining a new variable of integration y' = Z\j2dt and substituting in (7, we get 

r(^l9^ r + f 1 2s,„ + ,',--2Co(„ + „r-. 

+ 2I)a^{a + y'f-°^ -{2B -2C + 2D)a'-'^\} (16) 

It turns out that one has to be very careful to extract the asymptotic behaviour of the above integral, because the 
correlations at scales y' < a and y' > a contribute different power laws, and they compete with each other. It is best 
to split the integral in Eq. ^| into two parts and . Accordingly (7 = Ci + C2 . In Ci we expand the terms 
of the integrand in powers of y' /a (with y' < a), while in C2 we expand the terms in powers of a/y' (with y' > a). 

The Ci integral, with its integrand expanded in a series yields: 

= T{d/2)i2dr/^ X ^ ^ t'~T57 + ^(^)] 

where 5i — Y{dj^^2d)^^l ~ i5^rf+i) +''']■ obtaining Eq. 1171 we have used the fact that for small y' and very large 

i, we may approximate e^^ Z^*^* 1 inside the integral. In Eq. El the factor a'*^" smoothly goes to zero if a 
for a < d. If we wish, by taking smaller and smaller cutoff a, we can actually ignore the contribution of of Ci with 
respect to C2. On the other hand, if a > c?, a ^ will lead to a divergence, and hence we have to keep a finite, and 
thus the contribution of Ci cannot be ignored. 

Now we turn to the integral C2 and it goes as follows: 

/lj-2 — d/2 roo Ac)i Ai2 2 n^-a 

« ao t + ai i2-"/2 + 0(t3/2-"/2) ^ ^^^d-a ^2-d/2 (^8) 

where the newly defined constants are oq = 4A/[r(d/2)(d - 2)], ai = i2Bd'^ /[T{d/2){d - a){2d)°'/\ 62 = 
P^{-4/[(d- 2)(1- a)(2- a)] - 8BdV(d - a) + 0(a) + 64}, and 64 = 24/[(d- 4)(1 - a)(2 - a)(3- a)(4 - a)]. 
In going from the first line to the second line in Eq. I18l we have made the approximation of replacing the upper limit 
00 by \/2di, as the exponential factor exp{—y''^/2dt) sharply cuts off the integral beyond that scale anyway. After 

that we have set e~'^ Z^''* ~ 1 and then done the integrals. We see that in Eq. 1181 C2 has three dominant power laws 
with powers equal to 1, (1 — a/2), and (1 — d/2). Combining Eqs. 1171 and lT^ we get finally the two-time incremental 
correlation function from Eq. 1151 as 

C(ti,i2) -C(t) «Bii + ait2-"/2 + S2<2-'i/2, (19) 

where Bi = 1 -I- ag, and B2 = (&i -I- b2)a'^~°'- This implies that asymptotically, at long time t, C{t) ^ t^^ . Thus x{t) 
is a fBM process and the associated Hurst exponent H is given by 

i/ = max[l/2,(l-a/4),(l-d/4)] (20) 
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H= 1 - a/4 



H= 1/2 



H = 1 - d/4 
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FIG. 1; The values of the exponent H in the three regions of the parameter space of d versus a are shown. 



More explicitly, it means that for d < 2, for a < d, = (1 — a/4) while for a > d, H = {1 — d/A). On the other 
hand, for d> 2, for a < 2, H — {1 — a/A) while for a > 2, H — 1/2. This scenario is shown in Fig. ^ 

We have analysed the cases of marginal behaviour separately. For short-range disordered MdM model, d = 2 was 
shown \n\ to have marginal behaviour. For the current model, we have done explicit calculations (similar as above) 
for the cases a = 2, d = 2 and a = d [2l|. (i) For a = 2, for large t, C{t) ~ t In i (for d > 2), ~ ^^^^ 
d < 2), and ^ t (In i)^ (for d = 2). This means that the semi-infinite line a — 2 with d>2 (see Fig. has marginal 
behaviour, (ii) For d = 2, C{t) ~ Hn t (for a > 2) and - ^2-0/2 (■^q^. q, ^ 2). Hence the semi-infinite line d = 2 with 
a > 2 (Fig. ^ exhibits marginality. (iii) Finally we also find logarithmic correction, namely C{t) ~ t^~'^/^ln t, along 
the hne a = d, for a and d both < 2. 

From the known first-passage property of fBm mentioned earlier, the disorder and thermal averaged probability 
that the process x{t) does not cross zero up to time t decays as a power law, P{t) ^ t^^ with 6 — 1 ~ H. Using the 
results for H in Eq. (|20|l . we get 

61 = min[l/2,a/4,d/4]. (21) 

This means that for d < 2, if a < d, = a/i while if a > d, 6* = d/4. For d > 2, if a < 2, = a/A, while if a > 2, 
e = 1/2. 

Although for d < 2 and a > d, we expect analytically C{t) ~ ^2-d/2 large t, we find numerically (see section 
IIVI) that it is rather tricky to find this power law. The latter is a consequence of the Eq. 1171 which tells us that this 
long time behaviour arises, predominantly, as a contribution of correlations at distance scales y' < a. On the other 
hand Eq. 1181 shows that the power law ^2-0/2 g^j-jges out of a contribution due to correlations at scale y' > a. So if a 
finite cutoff a is not kept, and the distance scales y' <^ a are not probed carefully, numerically one may find a wrong 
asymptotic power namely 2H = (2 — a/2) instead of the expected power 2H = (2 — d/2). We will discuss this further 
in the next section Hvl 



IV. NUMERICAL SIMULATION AND RESULTS 



Since the exponents in Eq. 1201 followed from Eq. 1191 which is an approximate asymptotic form of the exact Eq. 1151 
we first checked by numerically integrating Eq. ^jthat the predicted powers are indeed true. While doing this, we 
realized that in order to see the power C{t) ~ ^2-^/2 ^^le regime a > d (for d < 2) one should be very careful to 
choose the numerically finite integration measure SZ <C 1 (where 6Z approximates infinitesimal dZ), in order to see 
the asymptotics. This is related to the comment made earlier that we have to probe y' <C a (where y' = Z\j2dt as 
mentioned earlier) to get the asymptotic contribution, in this regime. 

We now turn to verifying the analytical predictions of the previous section, namely H as in Eq. I2UI and the 
corresponding exponent Q as in Eq. 1211 by doing numerical simulations. We use the time discretised version of Eqs. 
Hand El 

yiiUn+i) = yt{t,n) + VAirii{t„i) (22) 
x{tm+i) = x{t„-,) + Atv {x{t,n)) + y^r]^{t^), (23) 
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FIG. 2: Log-log plot of C'{t) versus t for values of a = 0.6 (above) and 0.8 (below). The thin dotted straight lines are fits to 
the data (in thick lines) with a slope equal to 2H — 2 — a/2 as given by Eq. 1201 and at large times t the data approach the 
straight lines. The data curves are scaled by suitable factors for visual clarity. 
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FIG. 3: Log-log plot of C{t) versus t for values of a = 3.5 (above) and 1.5 (below). The thin dotted straight lines are fits to 
the data (in thick lines) with a common slope equal to 2H = 3/2, as in Eq. 1201 Both the curves approach this slope at large 
times. The curves are scaled for visual clarity. 



where tm = mAt, with m being an integer. We choose At < 0.5 in our simulations so that the stabihty is guaranteed 
|16|. The variables rii{tm) and rix{m) are independent Gaussian variables for all tm and each is distributed with zero 
mean and unit variance. Moreover our simulations here are all for d = 1, i.e. we have a single transverse coordinate 
y. Hence the subscript i assumes a single value equal to 1. 

Unlike pjlll5j|. here we had to numerically create quenched power law correlated velocity v{y) of the form given by 
Eq. O We note that if we have a delta-correlated velocity correlation in Fourier /c-space with a fc-dependent prefactor 
as follows 

i{ki)i{k2) = c{ka)\ki\°'-^S{ki + fe), (24) 

with constant c(fca) = [cos(fca) dp^^^^ + sin(A:a) dp^^] , in the real y-space we have v{yi)v{y2) given by Eq. 
IHup to an constant factor. The factor c{ka)\ki\"~^ in Eq. |^is the square of the variance of the distribution of v 
in fc-space. In our numerics, we first generate a realization of Gaussian distributed random numbers v{k) in fc-space 
with zero mean and variance proportional to ^/c{kaj\k\^°'~^'>^^ , and then Fourier transform back to real space. The 
latter realization of random v{y) are thus guaranteed to have the required power law correlation in real y-space. 

We note that to see the correct asymptotic behaviour, one has to explore distance scales y <C a. We consider a grid 
along the y-direction with grid spacing Ay, extending from —N = —ymax/ Ay = —500000 to = ymaxi Ay = 500000. 
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FIG. 4: Log-log plot of P{t) versus t for the values of a = 0.75 (above) and 3.5 (below). The straight lines are fits to the data 
with a slope equal to —9 where 9 = 1 — H (see Eq. 



At each point of this grid we choose v(jj) by a discrete Fourier transform of the random numbers v{km) obtained for 
all km (with km = —N, . . . , N). Once a set of {v{y)} is thus chosen, they remain fixed at all times during different 
thermal histories. This set {v{y)} constitutes a particular realization of disorder. Finally one performs the disorder 
average (. . .) by averaging over various realizations of the set {v{y)}. In all our simulations we have averaged over 
100 realizations of disorder using 1000 thermal realizations for each set of disorder. 

There are further, three important things to be noted. Firstly, for a < d, the role of cutoff a is unimportant and 
hence it may be safely ignored. We have verified for d = I that the results are unchanged in comparison to finite a, by 
taking a — > 0, which is equivalent to setting c{ka) = 1 in Eq. 1241 However, the role of a in the a > d regime cannot 
be overemphasized. The cutoff a has to be chosen such that a > Ay. We have used a = 10 Ay in our simulations. 
For low values of the ratio a/ Ay , the expected behaviour of C{t) ^ jg j^q^ seen. The second point is that the 

value k = leads to numerical difficulties for a < d. The latter was taken care of by using a very small finite value 
of k in place of A: = 0. Finally we note that the two integrals in the constant c{ka) are rapidly decaying integrals and 
better results are obtained by using very small grids of integration. 

In figure 2, we see that for different a < d{= 1), the incremental correlation function C(t) ^ t^^ with exponent 
H = 1 — a/4, as predicted by Eq. [5D| In figure 3, we see that for different a > d{= 1), again C{t) ^ t^^ , but now 
with 2H = 2-d/2 = 3/2. 

We directly measured the averaged persistence probability P{t) in our simulations. The procedure is described 
clearly in |1^. Our results are shown in figure 3. We have shown P{t) for two different a values: for a < d{~ 1), the 
power law exponent = a/4 as expected, while for a > d{= 1) we have 6 = d/A =1/4. 

V. REMARKS AND CONCLUSION 

In this paper, we have derived the incremental two-time correlation function C{t) and persistence probability P{t) 
of the longitudinal x-coordinate of a random walker subjected to a quenched random power law correlated velocity 
field. This problem is a variant of the original MdM model, now with long range velocity correlation. The main result 
is shown in figure ^ There are three distinct regions in the parameter space of d versus a. (i) For low dimensions 
{d < 2) and large power law exponents {a > d), the behaviour of C(t) and P{t) are same as that of the short-range 
correlated MdM model !l l1 . Thus this regime is short-range disorder dominated, (ii) On the other hand for strong 
long-range disorder i.e. a < min[(i, 2], the temporal behaviour of C{t) and P(t) depend on a continuously, so that 
long-range correlation really has a new effect. The latter regime is long-range disorder dominated. {Hi) Finally for 
d > 2 and a > 2, the Hurst exponent is 1/2 , i.e., the same as that of a simple random walker (SRW). Quenched 
disorder is irrelevant, and this regime is thus thermal noise dominated. Further, at the boundaries in the Fig. ^ 
denoted by the dotted lines, there are logarithmic corrections to the power law expressions for C{t) and P{t). 

We note that the the quenched random velocity field is having an effect on the temporal behaviour of a SRW, 
only when the disorder is sufficiently strong and transverse dimension d is small. For short-range correlated MdM 
model pj| the H changes only when the number of transverse directions d < 2. It is as if the velocity kicks in the 
x-direction fail to have much impact when the particle has many more transverse directions to explore. In this paper 
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we find that for a < min[c?, 2], i.e. for stronger correlated or longer ranged disorder, the short-range correlated MdM 
model exponents give way to new exponents. Thus even if transverse dimensions are large, sufficiently long-range 
disorder makes the random walker move super-diffusively. For high dimension {d > 2) and weak disorder (a > 2), the 
quenched noise looses out to the thermal noise, and the asymptotic temporal properties are that of simple Brownian 
motion. 

Wo also note that the short-distance cutoff scale a of the power law, plays a crucial role in this problem. For the 
case a > d, numerically, one has to be very careful to explore distance scales much smaller compared to the power-law 
cutoff scale a. It is one of the most important aspects of this regime that very short distance correlations seem to 
control the long time behaviour. For the regime oi a < d, since correlations at scale y' > a determine the asymptotic 
behaviour rather than y' < a, the cutoff a may be ignored. 
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